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ABSTRACT 

We present an analysis of Spiteer/IRAC primary transit and secondary eclipse 
lightcurves measured for HD 209458b, using Gaussian process models to marginalise 
over the intrapixel sensitivity variations in the 3.6^m and 4.5/im channels and the 
ramp effect in the 5.8/im and 8.0^m channels. The main advantage of this approach 
is that we can account for a broad range of degeneracies between the planet signal 
and systematics without actually having to specify a deterministic functional form 
for the latter. Our results do not conhrm a previous claim of water absorption in 
transmission. Instead, our results are more consistent with a featureless transmission 
spectrum, possibly due to a cloud deck obscuring molecular absorption bands. For the 
emission data, our values are not consistent with the thermal inversion in the day- 
side atmosphere that was originally inferred from these data. Instead, we agree with 
another re-analysis of these same data, which concluded a non-inverted atmosphere 
provides a better ht. We Hnd that a solar-abundance clear-atmosphere model without 
a thermal inversion underpredicts the measured emission in the 4.5^m channel, which 
may suggest the atmosphere is depleted in carbon monoxide. An acceptable ht to the 
emission data can be achieved by assuming that the planet radiates as an isothermal 
blackbody with a temperature of 1484 ± 18 K. 

Key words: planets and satellites: atmospheres; planets and satellites: general; meth¬ 
ods: data analysis, astronomical instrumentation, methods, and techniques; stars: in¬ 
dividual: HD209458 


1 INTRODUCTION 

Over the past decade, the Spitzer Space Telescope has proven 
to be a productive facility for characterising the atmospheres 


of transiting exoplanets (e.g. Charbonneau et al. 


2005 


Dem- 


ing et al. 112006 

[Knutson et al. 2007a 

Desert et al.| 

20091 

Crossfield et al. 

2012 

Lewis et al. 2013 

Todorov et al. 2014 

1. 


The ability of its instruments to probe the ~3-25/rm wave¬ 
length range has provided constraints on the thermal emis¬ 
sion of numerous exoplanets, as well as atmospheric trans¬ 
mission in a region dominated by absorption from molecular 
species such as water, methane, carbon monoxide, and car¬ 
bon dioxide. Hot Jupiters have offered especially favourable 
targets for such observations, given their large atmospheric 
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scale heights and relatively strong emission at these wave¬ 
lengths. 

This paper focuses on observations made with the In¬ 
frared Array Camera (IRAC), which has been the most 
widely used Spitzer instrument for observing exoplanets. 
Specifically, we analyse ten transits and eleven eclipses that 
have been measured for HD 209458b, most of which have al¬ 
ready been published (Knutson et al.||2008 Beaulieu et al. 


2010 Zellem et al. 2014 Diamond-Lowe et al. 20141. By 
providing a uniform analysis of these datasets, we aim to 
re-evaluate a number of claims that have been made in the 


literature. In particular, Beaulieu et al. (20101 measured sig¬ 


nificantly larger effective radii for the planet in the 5.8pm 
and 8.0pm channels relative to the 3.6pm and 4.5pm chan¬ 
nels in transmission, and interpreted this as evidence for 


water absorption. However, Deming et al. (20131 have since 


resolved the water absorption band centred at 1.4pm us- 
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2 Evans et al. 


ing the Hubble Space Telescope (HST) Wide Field Camera 
3 (WFC3), and found that it has a much lower amplitude 
than would be expected based on the results of Beaulieu et 
al. We also seek to address the claim of a thermal inver¬ 
sion in the dayside atmosphere, which was hrst postulated 


by Knutson et al. (20081 based on the deeper eclipses those 


authors measured for the 4.5/rm and 5.8/rm channels rel¬ 
ative to the 3.6/im channel. The former channels coincide 
with absorption features due to water and carbon monox¬ 
ide: therefore, seeing these features in emission would sug¬ 
gest an increasing temperature profile with decreasing pres¬ 
sure. Diamond-Lowe et al. (20141 have challenged this pic¬ 


ture by presenting revised eclipse depths that are sugges¬ 
tive of a non-inverted pressure-temperature profile. Further¬ 
more, Hansen et al. (20141 have suggested that the emission 


data for HD 209458b are consistent with radiation from an 
isothermal blackbody. 


The conflicting results obtained by different authors 
analysing the same datasets is likely due to the various 
methods that have been used to account for the instru¬ 
mental systematics that dominate IRAC lightcurves. The 
main contribution of the current study is to apply the ma¬ 
chinery of Gaussian processes (GPs) to the task of treating 
these systematics. This work follows similar applications of 
GP models to transit lightcurves published by [Gibson et al.| 
( 2012b|a 2013a|b I, Evans et al. (20131, and Gibson (20141. 
One of the primary advantages of GP models is that they al¬ 
low us to naturally handle correlations in the data that may 
be poorly understood from a hrst principles standpoint, by 
specifying only high-level properties of the covariance. This 
relaxes the assumptions built into our model, by removing 
the need to associate the systematics with a deterministic 
functional form. The resulting model is less restrictive, al¬ 
lowing us to capture a broad range of systematics behaviours 
with a relatively small number of tunable parameters. Fur¬ 
thermore, GPs are Bayesian models in the sense that uncer¬ 
tainty is treated transparently using self-consistent rules of 
probability. Each unknown in our model is associated with 
a probability distribution that rehects our uncertainty in 
its value, and it is possible to write down an expression for 
the likelihood of the observed data given specihc values for 
the model parameters, i.e. the model posterior distribution. 
We can then optimise the model posterior with respect to 
the unknown parameters, or marginalise over the parameter 
space using a method such as Markov chain Monte Carlo 
(MCMC). These properties make GP models suitable for 
inferring planet parameters from transit lightcurves affected 
by systematics that are not especially well-understood, such 
as those obtained with IRAC. 


The paper is arranged as follows. Section describes 
the lightcurve observations analysed for this study and Sec¬ 
tion describes how we produced lightcurves from the raw 
data frames. IRAC instrumental systematics are described 
in Section with an overview of methods that have been 
used to correct for them previously in the literature. Section 
[^outlines the GP methodology that we adopt in the current 
study, and Section [^describes the lightcurve fitting. The re¬ 
sults are presented in Section [T] and discussed in Section]^ 
with a focus on the implications for the planet atmosphere. 
Our conclusions are summarised in Section]^ 


2 OBSERVATIONS 


We have analysed ten primary transits and eleven secondary 
eclipses for HD209458b, made across all four IRAC chan¬ 
nels. Details of the relevant Spitzer observing programs are 
given in Table along with references to previously pub¬ 
lished analyses. More specihc information for the individual 
lightcurves is given in Table For this study, we did not 
model the complete half- and full-phase datasets acquired 
for Programs 30825, 40280, and 60021. Instead, for these 
lightcurves only ~5 hr subsections centred on the transits 
and eclipses were analysed. 

Most observations were made in stare mode. The only 
exceptions were those made as part of Program 20523, which 


have been published in Knutson et al. (20081. For the latter. 


four sets of 64 frames were acquired in a given channel be¬ 
fore the telescope was repointed to be centred on the next 
channel, and another four sets of 64 frames were taken. This 
process was cycled through each of the four channels, and 
repeated for the duration of the observations. Knutson et al. 
discarded the first set of 64 frames in each set of four, as the 
star was still drifting significantly during this time following 
the repointing. In addition, those authors discarded the first 
ten frames and the 58th frame from each set of 64 for the 
5.8/rm and 8.0pm channels, as these exhibited count levels 
consistently below the median. We performed two separate 
analyses for these lightcurves: one with this culling applied, 
and one without. However, we obtained consistent results in 
both cases, and only present results for the unculled dataset 
below. 


3 DATA REDUCTION 

The Basic Calibrated Data (BCD) frames for each lightcurve 
were reduced using a custom pipeline written in the Python 
programming language]^ The first step performed by the 
pipeline is to calculate the background level and locate the 
stellar centroid in each BCD frame. The background was 
estimated by taking the median pixel value from the four 
8x8 pixel subarrays at the corners of each frame, and then 
subtracted from each pixel in the array. The centroid co¬ 
ordinates were then determined by taking the flux-weighted 
mean of a 7 x 7 pixel subarray centred on the approximate lo¬ 
cation of the star. An initial guess was provided for the stel¬ 
lar centroid coordinates in the first frame, and coordinates 
determined for the previous frame were used as the initial 
guess in subsequent frames. Mid-times were computed for 
each exposure in Barycentric Julian Date Coordinated Uni¬ 
versal Time (BJDutc) using the BMJD_OBS and FRAM- 
TIME header entries. Bad frames were flagged by identifiy- 
ing those with outlying centroid coordinates or pixel counts. 
This was done by comparing against the median and stan¬ 
dard deviation of the 30 frames immediately preceding and 
following each frame. If the centroid coordinates or any pixel 
counts within a subarray spanning the photometric aperture 
centred at the stellar centroid differed from the median by 
> 5cr, the frame was discarded from the analysis. This pro¬ 
cess was iterated twice, resulting in 0.2-4.4% of the frames 
being discarded depending on the dataset (Table [^. 

^ Publicly available at www.github.com/tomevans 
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Table 1. HD209458 datasets analysed for this study. 


Channels (/im) 


Program 

P.I. 

Type 

20523 

D. Charbonneau 

Eclipse 

40280 

H. Knutson 

Half-phase 

461 

G. Tinetti 

Transit 

60021 

H. Knutson 

Full-phase 


3.6 

4.5 

5.8 

8.0 

Reference^ 

Yes 

Yes 

Yes 

Yes 

Kn08, DL14 

- 

- 

- 

Yes 

DL14 

Yes 

Yes 

Yes 

Yes 

BelO 

Yes 

Yes 

- 

- 

Zel4, DL14 


BelO ( Beaulieu et al.|2010 l, DL14 l Diamond-Lowe et al.|2014[ |, Kn08 ( Knutson et al.|2008[ l, Zel4 ( Zellem et al.|2014 l. 


Table 2. Lightcurve details. 


Program Type Channel Date 

Gm) (UT) 


Mode“ 


Flaggeq 

(%) 


3 ivj] 




(sec) 


(pijg 


20523 

Eclipse 

3.6 

2005 Nov 28 

sub, 0.1 

0.36 

1115 

32 

4.5 

4.0 


Eclipse 

4.5 

2005 Nov 28 

sub, 0.1 

0.21 

1115 

32 

4.5 

3.0 


Eclipse 

5.8 

2005 Nov 28 

sub, 0.1 

0.32 

1115 

32 

4.5 

2.5 


Eclipse 

8.0 

2005 Nov 28 

sub, 0.1 

0.55 

1115 

32 

4.5 

3.0 

40280 

Transit 

8.0 

2007 Dec 25 

sub, 0.4 

1.11 

1577 

32 

13.5 

3.5 


Eclipse 

8.0 

2007 Dec 24 

sub, 0.4 

1.12 

1577 

32 

13.5 

4.0 

461 

Transit 

3.6 

2007 Dec 31 

full, 0.4 

1.68 

1427 

2 

16.8 

2.5 


Transit 

3.6 

2008 Jul 19 

full, 0.4 

2.03 

1428 

2 

16.8 

2.5 


Transit 

4.5 

2008 Jul 22 

full, 0.4 

1.67 

1288 

2 

16.8 

2.5 


Transit 

5.8 

2007 Dec 31 

full, 2.0 

3.78 

1425 

2 

16.8 

3.5 


Transit 

5.8 

2008 Jul 19 

full, 2.0 

3.81 

1425 

2 

16.8 

3.0 


Transit 

8.0 

2008 Jul 22 

full, 2.0 

4.42 

1286 

2 

16.8 

4.0 

60021 

Eclipse (1st) 

3.6 

2011 Jan 12 

sub, 0.1 

0.25 

1285 

128 

16.8 

2.5 


Transit 

3.6 

2011 Jan 14 

sub, 0.1 

0.25 

1286 

128 

16.8 

2.5 


Eclipse (2nd) 

3.6 

2011 Jan 16 

sub, 0.1 

0.29 

1286 

128 

16.8 

3.0 


Eclipse (1st) 

3.6 

2014 Feb 13 

sub, 0.1 

0.23 

1286 

128 

16.8 

3.0 


Transit 

3.6 

2014 Feb 15 

sub, 0.1 

0.23 

1286 

128 

16.8 

3.0 


Eclipse (2nd) 

3.6 

2014 Feb 17 

sub, 0.1 

0.18 

1285 

128 

16.8 

2.5 


Eclipse (1st) 

4.5 

2010 Jan 18 

sub, 0.4 

0.62 

1577 

32 

13.6 

2.5 


Transit 

4.5 

2010 Jan 19 

sub, 0.4 

0.73 

1575 

32 

13.6 

3.0 


Eclipse (2nd) 

4.5 

2010 Jan 21 

sub, 0.4 

0.62 

1577 

32 

13.6 

2.5 


“ Readout mode and frame time in seconds. 

^ Fraction of frames flagged as bad. 

^ Number of frames in the binned dataset used for the GP lightcurve fit. 

Number of consecutive frames per bin. 

® Median cadence of binned frames. 

/ Photometric aperture radius. 


Photometry was performed for each remaining frame 
by summing the pixel counts within circular apertures. Sep¬ 
arate reductions were obtained for different aperture radii, 
ranging between 2-6 pixel in increments of 0.5 pixel. Due to 
the undersampled nature of the IRAC point spread func¬ 


tion (PSF), we linearly interpolated the native pixel array 
onto a 10 X 10 super-sampled grid, as has been done by 
others previously (e.g. Stevenson et al.||2010 |. These inter¬ 
polated sub-pixels were counted towards the aperture sum 
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if their centres fell within the aperture radius. The resulting 
lightcurves are shown in Figure 


works by linearly decorrelating the measured fluxes against 
the individual pixel counts within a subarray centred on the 
PSF, rather than using the centroid coordinates directly. 


4 INSTRUMENTAL SYSTEMATICS 

The raw lightcurves are affected by instrumental systematics 
that are characteristic of IRAC and have been documented 
extensively in the literature (e.g. Charbonneau et al.||2005| 
Ago! et aL]|2010| [Seager fc Deming||2010| [Stevenson et al'| 
2012). The systematics divide into two broad categories: in¬ 


trapixel sensitivity variations in the 3.6/im and 4.5/im chan¬ 
nels, and the ramp effect in the 5.8pm and 8.0pm channels. 


4.1 3.6pm and 4.5pm channels 

In the 3.6pm and 4.5pm channels, which employ InSb de¬ 
tectors, the measured flux correlates with the position of the 
stellar PSF on the detector array. As noted in the IRAC in¬ 
strument handbook|^this effect is believed to be caused by 
variations in the quantum efficiency across individual pixels. 
Pointing drift during observations, combined with the un¬ 
dersampled PSF, therefore results in variations in the mea¬ 
sured flux at the level of a few percent. 

Traditionally, intrapixel sensitivity variations have been 
treated in IRAC data by decorrelating the lightcurve against 
a low-order polynomial in the xy centroid coordinates of 
the stellar PSF, which can either be removed prior to fit¬ 
ting the planet signal or ht simultaneously with the planet 
signal (e.g. [Charbonneau et al.||2008| [Knutson et al.|[2008[ 
Desert et al. 2009). The main issue with this method is 


4.2 5.8pm and 8.0pm channels 

In the 5.8pm and 8.0pm channels, which employ Si:As de¬ 
tectors, the measured flux smoothly increases or decreases 
before levelling off. As with the intrapixel sensitivity varia¬ 
tions, the amplitude of this effect is usually a few percent, 
with the steepest change in flux occurring during the initial 
~ 1 hr of observations. The ramp can be attributed to the 
sensitivity of individual pixels varying as a function of time, 
where the rate of change depends on the illumination level 
of the pixel (e.g. Knutson et al.|2007a l. 

It has been suggested that the ramp is caused by elec¬ 
trons getting caught in charge traps within the detector (e.g. 


Deming et al. 2007 Agol et al. 2010 Seager & Deming 20101. 


As more photons arrive, the charge traps hll up, resulting in 
less electrons getting trapped and a higher flux being mea¬ 


sured. However, Seager & Deming (2010) acknowledge that 


that a low-order polynomial may not have the flexibility to 
fully capture the underlying correlations in the data if there 
is fine-scale structure present. This could be addressed to 


the physics of the detectors are poorly understood and this 
simple picture may be incomplete or wrong. For instance, 
lightcurves measured in the 5.8pm channel exhibit a ramp- 
down behaviour (Figure]^, which is not obviously explained 
by charge trapping. 

Standard practise has been to fit 5.8pm and 8.0pm 
lightcurves by multiplying the transit signal by a paramet¬ 
ric model that provides a good approximation to the ramp. 
One approach is to remove the initial steep section of the 
lightcurve and fit the remainder with a linear or quadratic 
polynomial in time (e.g. [Deming et al.|2007||Beaulieu et al.j 
2010). However, a more common approach is to model the 


full baseline with a low-order polynomial in logarithmic time 


some extent by continuing to add higher order terms to the 

(e.g. |Charbonneau et aL||2008| |Knutson et aL||2008| |Desert| 

polynomial decorrelation; however, such an approach runs 

et al.[|2009 

Machalek et al.[[2010 1 or exponential time (e.g. 

the risk of overfitting and increases the dimensionality of 

Agol et al. 

20101. 

the parameter space that must be marginalised over using a 

Subtle differences between the various ramp functions 


method such as MCMC. 

An alternative decorrelation method has been suggested 


by Ballard et al. (20101, which uses a 2D Gaussian convolu¬ 


tion to construct a smoothed pixel sensitivity map from the 
flux measurements. Using this approach, the authors identi¬ 
fied a high-frequency corrugation structure in the xy sensi¬ 
tivity for a 4.5pm lightcurve that would not have been cap¬ 
tured by a low-order polynomial decorrelation. Along similar 


lines, Stevenson et al. (2012) have proposed the use of bilin¬ 


ear interpolation of the measured fluxes to map intrapixel 
sensitivity variations, which is sensitive to spatial scales cor¬ 


responding to the grid of knots used. Knutson et al. (20121 
and [Lewis et al. ] ( [Ms] ) have also identified a correlation 
between the noise pixel parameter /3 of a given data frame 
- which is inversely proportional to the PSF sharpness pa¬ 
rameter, described in [Muller fc Buffington| ( [l974[ ) - and the 
measured flux. This can be treated by either using vari¬ 
able photometric apertures for each frame that scale with 
P, or including /3 as a decorrelation variable in the system¬ 


atics treatment. More recently, Deming et al. (2014) have 


presented a pixel-level decorrelation (PLD) method, which 


http://irsa.ipac.caltech.edu/data/SPITZER/docs/irac 


used to model 5.8pm and 8.0pm lightcurves can potentially 
bias the inferred planet parameters. For example, polynomi¬ 
als in exponential time tend to retain slightly more curvature 
than polynomials in logarithmic time after the initial steep 
gradient. This can result in an underestimated transit depth 
if the ramp is decreasing, or vice versa if the ramp is increas¬ 
ing. To try avoid effects such as these, the lightcurve analysis 
is often performed separately for a number of different ramp 
functions, and the one that minimises the model residuals 
or maximises an approximation to the Bayesian model ev¬ 
idence, such as the negative Bayesian information criterion 
([Schwarz|1978|), is retained (e.g. [Stevenson et al.[[201^. 


5 GAUSSIAN PROCESSES FOR IRAC 
SYSTEMATICS 

In this paper, we present an alternative method for treating 
intrapixel sensitivity variations and the ramp effect in IRAC 
lightcurves based on Gaussian processes (GPs). The first ap¬ 
plication of GP models to transit lightcurves was made by 


Gibson et al. (2012b I, with subsequent applications in Gib 
son et al. (2012a 2013a|b ), Evans et al. (2013), and Gibson 
(2014), to which the reader is referred for further details. 
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Transits 


(a) 3.6/;m 2007 Dec 31 


1.15 


1.10 


1.05 


XI 

a 

QJ 

> 


1.00 


0.95 


0.90 


0.85 



-4-3-2-10 1 2 3 

Time from mid-transit (hrs) 


1.15 


Eclipses 


1.10 


1.05 


1.00 


0.95 


0.90 


0.85 


(a) 3.6/im 2005 Nov 28 





(b) 3.6/im 2011 Jan 12 


(c) 3.6/(m 2011 Jan 


(d)3.6//m 2014 Feb 13 


(e)3.6//m 2014 Feb 17 



(f) 4.5/nn 2005 Nov 28 



(h) 4.5/m 2010 Jan 21 




(i) 5.8/m 2005 Nov 28 , 






(k) 8.0/m 2007 Dec 25 


-4-3-2-10 1 2 3 4 

Time from mid-eclipse (hrs) 


Figure 1. Raw lightcurves obtained for the 10 transits (left column) and 11 eclipses (right column) analysed in this study. Red lines 
show best-fit GP models, which are described in Section]^ The lightcurves shown in this plot have been binned in time, which was done 
to make the GP model fitting computationally tractable (Section |5.3[ l. The plotted GP models have been further binned into 1 minute 
bins and vertical offsets have been applied to each lightcurve for clarity. 
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Formally, a GP is defined as a collection of data points, 
any subset of which has a multivariate normal distribution 
(e.g. Rasmussen fc Williams|2006 1 . Indeed, this assumption 
is regularly made for transit lightcurve analyses published in 
the literature, even if it is not stated explicitly. It is equiva¬ 
lent to assuming that what we measure is some underlying 
signal, which is the combination of the astrophysical signal 
of interest and additional systematic terms, plus uncorre¬ 
lated Gaussian noise. If we are in possession of a model /i 
that describes the underlying signal, the probability of mea¬ 
suring a specific dataset d = {di, d 2 , ..., dw} with Gaussian 
errorbars for a set of model 

parameters a = {qi, 012 ,..., ajv/} is given by: 


p(d|Q:) 


n 



(di - iiif 


= , 


( 1 ) 


where N denotes a multivariate normal distribution and 
S = diag[cri„] is the diagonal covariance matrix. For fixed 
errorbars, optimising the log likelihood is therefore identical 
to the familiar practice of minimising the statistic, since: 


11 "^ 1 

lnAf(/x,E) = -N\n2^ . (2) 

i=l 

Note, however, that Equationj^requires us to build all of our 
information about the underlying signal into the determinis¬ 
tic mean function /x. This is not desirable if, as is commonly 
the case, the systematics are poorly understood from first 
principles and an explicit functional form is not available to 
describe them. 

An alternative option is to incorporate the systematics 
into our model by allowing for nonzero off-diagonal entries in 
the covariance matrix of the likelihood function, such that: 


lnp(d|Q:, 7 ) = In Af (/x, K-|-E) , (3) 


where Kij gives the covariance between the ith and jth data 
points, and 7 are the parameters that control the behaviour 
of the covariance. Although we no longer have to provide 
an explicit functional form for the systematics contribution, 
we must now specifiy a kernel function to populate the en¬ 
tries of the covariance matrix K. However, by modelling 
the statistical covariance between data points rather than 
the deterministic systematics signal directly, it is possible to 
capture a broad range of behaviours with relatively few free 
parameters. Thus, GP models are simultaneously parsimo¬ 
nious and flexible. 

Before proceeding to describe the specific mean func¬ 
tions and covariance kernels adopted in the current study in 
Sections |5.1| and |5.2| it is worth also pointing out the fact 
that GPs have the desirable property of automatically im¬ 
plementing the principle of Occam’s razor. This can be seen 
if we expand Equationinto its constituent terms: 

lnp(djQ:, 7 ) = (K + S)"^ r - iln|K + E| - ]^N\n2'K , 


Table 3. Nonlinear limb darkening coefficients. 


Channels (/xm) 



3.6 

4.5 

5.8 

8.0 

Cl 

0.5564 

0.4614 

0.4531 

0.4354 

C2 

-0.5462 

-0.4277 

-0.5119 

-0.6067 

C3 

0.4315 

0.3362 

0.4335 

0.5421 

C4 

-0.1368 

-0.1074 

-0.1431 

-0.1816 


plexity of a model is equivalent to assigning similar probabil¬ 
ities to an increasing diversity of functions. In other words, 
as the model complexity increases, the likelihood function 
becomes less sharply peaked near the mean /x and the prob¬ 
ability mass of the model becomes more diffusely spread 
throughout the function space. This is precisely what hap¬ 
pens as the term — | ln|K -1- S| decreases, in effect penal¬ 
ising model complexity. Finally, the third term, — |A^ln27r, 
remains constant for a given dataset. The balance between 
the first two terms of Equation therefore ensures that the 
probability mass of the model is distributed over the param¬ 
eter space in a manner that optimises the trade-off between 
goodness-of-fit and model complexity. 


5.1 Mean functions 


The mean function /x defines the model for the astrophysical 
signal with well-understood form; namely, a primary transit 
or secondary eclipse. Eor this purpose, we adopt the analytic 


transit functions of Mandel & Ago! (2002 1 . We set the orbital 


eccentricity to zero, consistent with observational evidence 


(e.g. Pont et al. 20111. We fix the orbital period to P = 
3.52474859day ( Knutson et al.|2007b l. 

For the primary transits, the mean function parame¬ 
ters that are allowed to vary are the radius ratio Pp/P*, 
normalised semimajor axis a/P*, impact parameter b = 
a cos i/p*, and transit mid-time Tmid, such that 0 . = 
{Pp/P*,a/p*, 6 ,Tniid} in Equation]^ Stellar limb darken¬ 
ing is treated using the nonlinear law o f|Claret (2004 1 , with 
coefficients fixed to those provided by |Hayek et al. (20121 
which were obtained specifically for the IRAC bandpasses 
using a 3D stellar model for HD 209458 (Table |^. 

For the secondary eclipses, only the eclipse depth Pp/P* 
and eclipse mid-time Tmid are allowed to vary, such that ol = 
{Pp/P*, Tmid}- The remaining mean function parameters are 
fixed to the values published by Torres et al. {20081, namely. 


Pp/P* = 0.121, a/p* = 8.76 and b = 0.507. 


( 4 ). 2 


Covariance kernels 


where r is a vector containing the model residuals, with ith 
term given by Vi = di — fii. The first term on the righthand 
side, — |r^(K -f E)“^r, serves as a goodness-of-fit term. For 
a given covariance matrix, it increases as the residuals be¬ 
come smaller, rewarding mean functions /x that match the 
data well. The second term, — ^ In |K-hS|, can be thought of 
as a complexity penalty. This is because increasing the com- 


Entries of the covariance matrix K are constructed using a 
kernel function, such that Kij = fe(vi,Vj), where vi and vj 
are vectors of inputs associated with the ith and jth data 
points, respectively. By inputs, we refer to variables that 
correlate with the measured signal - these are typically the 
same variables that would be used for standard polynomial 
systematics decorrelations. Further discussion of common 
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GP kernels, such as the squared exponential and Matern 
kernels, can be found in |Gibson et al. 12012b I. 

In this study, we parameterise the entries of the covari¬ 
ance matrix K for the 3.6/im and 4.5/im channels as the 
sum of two kernels: a squared exponential kernel with the 
centroid xy coordinates as inputs and a Matern = 3/2 
kernel with time t as the input. Writing this out explicitly, 
the combined kernel is given by: 


k{vi, Vj) = k^y + kt , 

where: 


kxy — -^xy 6Xp 


Xi — Xj 
Lx 


kt — At 


1 + 


V3 


exp 


Vi - Vj 

Loj 


ti — tn 


V3 


(5) 


( 6 ) 


(7) 


such that 7 = { A^y , , Ly , At, Lt} in the notation of 

Equation]^ The squared exponential component of this ker¬ 
nel accounts for the smooth spatial variations in pixel sen¬ 
sitivities that dominate the systematics of the 3.6 /rm and 
4.5 ^m channels, while the Matern component accounts for 
any residual correlated noise in the lightcurve. 

For the 5.8/rm and 8.0/rm channel lightcurves, we used 
a squared exponential kernel to model the dominant ramp 
effect, with form given by: 


fcr = A^ exp 


( 8 ) 


where r = ln{t + h) is logarithmic time t, and /i is a pa¬ 
rameter that can be inferred from the data (see below). Pa- 
rameterising the covariance with k-r allows us to capture 
the dominant behaviour of the ramp effect; namely, a steep 
initial gradient followed by a levelling off of the measured 
flux (Figure [^. However, we stress that by parameterising 
the covariance between data points according to Equation 
we are not constraining the systematics to be monoton- 
ically increasing or decreasing in time. This is a valuable 
property of GPs given that the ramp effect is not necessar¬ 
ily strictly monotonic, with an overshoot effect identified in 
a number of datasets (e.g. Knutson et al.|2012 i. As with the 
other channels, we also include a time-dependent Matern 
1 / = 3/2 kernel kt to account for residual correlations in the 
lightcurve that may not be related to the ramp. Therefore, 
the final kernel is given by: 


k{vi, Vj) = k^ + kt 


(9) 


with covariance parameters 7 = { At , Lt , h j At , Lt} ■ 

The covariance kernels outlined above allow for system¬ 
atics treatments that are at least as versatile as others used 
in the literature. For instance, the k^y component of Equa¬ 
tion is similar in concept to the 2D Gaussian correction 
developed by Ballard et al. (Section 4.11. Similarly, the fer 


component of Equation is reminiscent of the logarithmic 
time polynomials used in other published studies (Section 
4.21. However, as it is only the covariance between data 
points that is parameterised in terms of logarithmic time, 
the underlying signal itself need not be monotonically in¬ 
creasing or decreasing. Furthermore, by parameterising the 
covariance rather than systematics signal directly, the GP 
model is capable of marginalising over a broader range of 
function space with relatively few tunable parameters. 


Before proceeding, we highlight the fact that the covari¬ 
ance kernels given by Equations [5]^ assume that the input 
variables are noise-free. While this is a reasonable assump¬ 
tion for time t, it is not necessarily the case for the centroid 
coordinates x and y. Indeed, the undersampled nature of 
the IRAC PSF makes the centroid estimates particularly 
susceptible to shot noise of individual pixels. However, due 
to the fact that this noise is white, and because in practise a 
small patch of a single pixel is densely sampled by the PSF 
over the course of a few hours, we expect its effect to be 
averaged out. For this reason, and in line with other anal¬ 
yses of IRAC lightcurves, we do not explicitly account for 
the noise of x and y. In future work, however, it may be 
worth considering a more explicit treatment of noisy inputs 


(e.g. Goldberg et al.|[l998 Mchutchon fc Rasmussen||2011 1, 
or simply smoothing noisy inputs before feeding them to the 
GP (e.g. Gibson et al.|2012b K 


5.3 Lightcurve binning 

GP models become computationally intractable for datasets 
with N > 10® data points. This is due to the need to fac¬ 
torise the N X N covariance matrix when evaluating the 
—/r®^(K + S)“®r term and computing the determinant 
— I In |K -f- S| for each log likelihood evaluation (Equation 
1^. Our code implements this using Cholesky factorisation, 
which has a computational cost scaling as O (). To apply 
GPs to IRAC datasets, most of which consist of > 10'* 
data points, we bin the fluxes and centroid xy coordinates in 
time prior to fitting. Binning factors were chosen according 
to the format of individual lightcurves, such that the time in¬ 
terval between successive binned points was < 15 sec and the 
number of binned points per lightcurve was Nb = 1000-1500. 
Lightcurves obtained in full array mode for Program 461 
were binned by a factor of two, giving a median cadence 
of about 17 sec; lightcurves obtained in subarray mode for 
Programs 40280 and 60021 with frame times of 0.4 sec were 
binned by a factor of 32, giving a median cadence of about 
14 sec; lightcurves obtained in subarray mode for Program 
60021 with frame times of 0.1 sec were binned by a factor 
of 128, giving a median cadence of about 17 sec; lightcurves 
obtained in subarray mode for Program 20523 with frame 
times of 0.1 sec were only binned by a factor of 32, due to 
the sparser sampling of the lightcurve as the telescope was 
constantly repointed during the observations, resulting in a 
median cadence of about 5 sec. Details are given in Table 
The obvious drawback of binning the lightcurves in time 
is that we lose information on the timescales of our bin sizes. 
For instance, the xy centroid coordinates can vary coher¬ 
ently over timescales < 15 sec (e.g. [Stevenson et al.|[20T^ . 
However, this should not affect our results significantly, as 
none of the astrophysical quantities of interest vary over 
the bin timescales, nor does the information content of the 
transit lightcurve degrade significantly as we reduce the time 
resolution to < 15 sec. Furthermore, high-frequency system¬ 
atics should mostly average out given the large number of 
binned data points for each tunable model parameter. Any 
correlations that remain will be accounted for by the time- 
dependent Matern ;/ = 3/2 kernel kt in our model (Equa¬ 
tions and . 
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6 LIGHTCURVE FITTING 


Lightcurves were fit individually by marginalising over 
the model posterior distribution using Markov chain 
Monte Carlo (MCMC), in order to quantify the degen¬ 
eracies between the astrophysical parameters of interest 
and instrumental systematics. Following Bayes theorem, 
the model posterior distribution is given by p(a, 7 |d) oc 
p(d|Q;, 7 )p(a)p( 7 ), where p(djQ:, 7 ) is the GP likelihood 
given by Equation and p{a) and p{'y) are the priors on 
the mean function and covariance parameters, respectively. 
We adopted uniform priors for the mean function param¬ 
eters a and covariance length scales Li. For the covari¬ 
ance amplitudes Ai, Gamma distribution priors of the form 
p{Ai) = Gam(l,100) oc exp[—lOOAi] for i = {t,xy,T} were 
adopted. The latter give decreasing probability to increas¬ 
ing covariance amplitudes, encouraging the GP to reduce 
the covariance amplitude unless justified by the data. The 
white noise level was included as a free parameter in 
each model, with a uniform prior. The ability to inflate the 
statistical errorbars above the formal shot noise floor pro¬ 
vides the models with some additional flexibility for dealing 
with high-frequency noise that may be present in the data, 
without having to reduce the correlation length scales Li of 
the covariance kernels to unreasonably small values. 


With the posterior distributions defined, the model fit¬ 
ting for each lightcurve proceeds as follows. Values for the 
model parameters were drawn randomly from the model 
prior, i.e. p{a)p{'y). With this as a starting point. Equa¬ 
tion [3 was optimised using the Nelder-Mead simplex algo¬ 
rithm ( [Nelder fc Mead|1965 1 to obtain maximum likelihood 
estimates (MLEs) for the parameters. A short Metropolis- 


Hastings MGMG chain (Metropolis et al. 1953 Hastings 


1970) of 1000 steps was initiated at the MLE, with step 


sizes pretuned to give acceptance rates of 20-40 %. The me¬ 
dian chain values were then used as the starting location 
for a second MLE optimisation. In practise, the randomness 
introduced by the short MCMC chains helped prevent the 
MLE optimisations from getting trapped at local maxima of 
the likelihood surface, thus increasing the chance of locat¬ 
ing the global likelihood maximum. To further increase this 
probability, the entire process was repeated ten times, each 
time from a different random starting point. This was done 
separately for each lightcurve produced using the different 
photometric aperture sizes (Section]^. The photometric re¬ 
duction giving the lowest scatter in the residuals was then 
selected for the remaining analysis. 


Before commencing the final MCMC chains, the covari¬ 
ance parameters were fixed to their MLE values, which are 
reported in Table [AT| This approach - which is often referred 
to as “type-II maximum likelihood” (for further discussion 


Gibson et al. 


2012 b[) - allows the expensive O ( 


covariance matrix factorisation required for the GP likeli¬ 
hood evaluation (Equation to be performed only once 
at the beginning of the chain. Subsequent steps only cost 
O ( N'^ ), resulting in much faster computations. The disad¬ 
vantage is that by Hxing the covariance parameters 7 , they 
are not marginalised over. In effect, this imposes an artifi¬ 
cial restriction on the range of systematics functions that 
are explored by the GP model. Consequently, there may 
be degeneracies between the planet signal and systematics 
that are not fully incorporated into the final uncertainties 


for the planet parameters a presented here. For example, 
in their re-analysis of the NICMOS transmission spectrum 
for HD 189733b, Gibson et al. (2012b 1 found uncertainties 
that were up to ~ 1.5 larger when the covariance parame¬ 
ters were allowed to vary in the marginalisation compared 
to when they were fixed to their MLE values. Therefore, the 
uncertainties we report in this study should be considered 
lower limits to the true uncertainties. 

It should be emphasised, however, that fixing the covari¬ 
ance parameters is quite different to hxing the parameters 
of an explicit functional model for the systematics. Instead, 
hxing the covariance parameters is somewhat analogous to 
selecting a family of parametric models for the systemat¬ 
ics, as they control the high-level properties of the function 
space spanned by the GP model. Rather than selecting from 
a handful of distinct parametric models, the GP model of¬ 
fers access to a continuum of possible functions. By using 
covariance parameters that optimise the GP likelihood, this 
continuum is narrowed in a principled manner. To compare 
with the bilinear interpolation method for pixel mapping 
used by Stevenson et al. (2012) for instance, optimising the 
covariance length scales L^ and Ly is similar to choosing the 
optimal grid spacing for the interpolation knots. Treating 
our dataset as a GP, we have the advantage of being able 
to do this in the context of a self-consistent probabilistic 
model, by maximising the likelihood function with respect 
to the unknown parameters using a numerical optimiser. 

Having fixed the covariance parameters, an initial chain 
of 10 ® steps was run with the planet parameters allowed 
to vary and step sizes again pretuned to ensure acceptance 
rates of 20-40 %. The first 5 x 10'* steps were discarded as 
burn-in. An additional four chains were then run for 10® 
steps each, with starting parameter values drawn randomly 
from normal distributions centred on the mean values of the 
hrst chain. The width of the normal distributions were taken 
to be five times the standard deviation of the first chain, to 
ensure the starting locations were well-dispersed in param¬ 
eter space. After discarding the hrst 5 x lO'^ burn-in steps 
of these chains, the Gelman-Rubin statistics for each pa¬ 
rameter were calculated (Gelman & Rubin 1992[ ). Eor all 
lightcurves, these were found to be well within 1 % of unity, 
consistent with the chains having reached stable states. Fi¬ 
nally, the hve independent chains were combined into a sin¬ 
gle chain, giving 2.5 x 10® samples from the posterior distri¬ 
bution. 


7 RESULTS 

Results of the primary transit and secondary eclipse MGMG 
analyses are given in Table Best-ht models are overplot¬ 
ted on the raw lightcurves in Figure and the corrected 
lightcurves with residuals are shown in Figure 


T.l Transmission (i?p/R*) 

We hnd good agreement in the inferred parameters across 
epochs for the majority of the transit lightcurves. The only 
exceptions to this are the 2011 Jan 14 and 2014 Feb 15 tran¬ 
sits measured in the 3.6/rm channel. For the 2011 lightcurve 
we obtain a value for Rp/Ri, that is > 6 cr discrepant relative 
to those obtained for the 2007 Dec 31 and 2008 Jul 19 3.6/im 
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Transits 




Time from mid-transit (hrs) Time from mid-transit (hrs) 
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Ph 
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Time from mid-eclipse (hrs) 


Eclipses 



(g) 4.5/im 2010 Jan 18 

* •. • •• 



(i)5.8/im2005 Nov 28 


(j) 8.0/ini 2005 Nov 28 
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4 
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a 


p. 

Ph 


a: 

13 

'a: 



-4-3-2-10 1 2 3 4 
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Figure 2. Corrected primary transit (top) and secondary eclipse (bottom) lightcurves obtained by dividing the raw lightcurves by 
the systematics components of the best-fit GP models. Model residuals are shown in the rightmost column, with grey shaded regions 
indicating the inferred white noise values a-w- 
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10 Evans et al. 


Table 4. Results of MCMC primary transit lightcurve analyses. Quoted values are the chain medians, with uncertainties corresponding to 
the ranges either side of the medians that contain 34% of the chain samples. Orbital inclination values i are derived from the corresponding 
impact parameter b = a cos i/R^ and normalised semimajor axis a/R~k values. Brightness temperatures Ti, are derived from the measured 
eclipse depths Fp /F* assuming an ATLAS stellar model | |Kurucz |1979[[T993[| for HD 209458 and integrating over the IRAC bandpasses. 


Transits 

Channel 

Date 

Tmid 

Rp/Rt, 

a/Ri, 

b 

i 

(pm) 


(BJDutc - 2450000) 




(deg) 

3.6 

2007 Dec 31 

4465.6371lto°°°27 

n 1 9n'77’”l"0-00085 
U.iZUf /_o.00084 

Q '7f)+0.26 

' ^-0.23 

n 597+0 035 
U.OZ/_o_o44 

86.541°;°° 


2008 Jul 19 

4666.54742t0;°°02i 

n 1 

U.izzzu_o,oo062 

o qq+0.33 

o.oy_0_26 

n 489+0 044 
U.452_q Qg2 

86.891°;®° 


2011 Jan 14 

5575.93102tg°°030 

U.1100^_Q QQQg7 

Q '7'7+0.33 

o- ' ' -0.33 

n t:9c:;+0.048 

U.OZO_Q Q5g 

86.571°;?® 


2014 Feb 15 

6703.85250t°;°°°ll 

u.iiyiy_o 00032 

Q 1 q + 0.10 

o.lci-o.io 

0.5901°;°?° 

85.841°;?? 

4.5 

2008 Jul 22 

4670.07290tO;°°029 

u.iziyy_Q QQQgj^ 

9 31+0-35 

y-J-L_o.30 

0.4231°;°®^ 

87.39l?|;?^ 


2010 Jan 19 

5216A0564tlZ°ool 

u.izuyy_Q QQQ29 

8.89tg;g? 

u.^yo_o 010 

86.821?;;?® 

5.8 

2007 Dec 31 

4465.63663t0;‘;°°i? 

U.iZUU/_o 00265 

9.20t0;|° 

0 435+0 066 
u.^oj_o.o82 

87.291°;®° 


2008 Jul 19 

4666.54744tg°°033 

n 1 iQQn+000284 

U.iiO»U_Q QQ272 

Q 9'7+0.38 

-0.37 

0 427+*^ *^^^ 

87.361?;;°? 

8.0 

2007 Dec 24 

4458.58730tg;°°g« 

'^•-‘-^^^^-0.00114 

8.78«;11 

n 51 3+0017 
U.0i9_o 019 

86.651?;;?? 


2008 Jul 22 

4670.07243t°;°°°i 

n 11991 

U.iiyyi_Q QQQyg 

8.67t0;?0 

n c;97+0 027 
U.oz/_o 032 

86.52l?;-|® 

Eclipses 

Channel 

Date 

Tmid 

Fp/F. 

n 



(pm) 


(BJDutc - 2450000) 

(%) 

(K) 



3.6 

2005 Nov 28 

3702.5274ltO;‘;°l^? 

u.uyci_o.o34 

144711°° 




2011 Jan 12 

5574.16955t°;°°;^g 

u.izz_o oil 

15911S 




2011 Jan 16 

5577.69695t°;°“04 

IJ--L^4_0.014 

1599l°® 




2014 Feb 13 

6702.09332t2;S 

^•-*^^^-0.017 

15451^° 




2014 Feb 17 

6705.61556tO;‘;°«6 

U.iUO_0.008 

1515l|? 



4.5 

2005 Nov 28 

3702.52942^0,0^1? 

n 145+0.018 

U.i40_o.oi7 

14741^2 




2010 Jan 18 

5214.64693^0;°°°®® 

U.i^U_Q QJ4 

145511^ 




2010 Jan 21 

5218.16915t0;0°;^^ 

n 1 ‘I'l+ooii 
U.iOO_0.011 

142611^ 



5.8 

2005 Nov 28 

3702.53033t0;00«7 

0 142+0'059 
^•-L^^_0.058 

1297li® 



8.0 

2005 Nov 28 

3702.5346lt0-0°°?2 

n 995+0.064 

U.zzo_o 063 

1433+219 

i^oo_2i6 




2007 Dec 25 

4460.35170t°;00?°° 

n 21 5+0012 
U.zio_o.oi2 

1397111 




lightcurves. This lightcurve has been classified as a failed ob¬ 
servation by the Spitzer Science Center due to the presence 
of high-frequency noise of unknown origin during the sec¬ 
ond half of the transit. For the 2014 lightcurve, we obtain 
values for a/R-^ and b that are both ~ 2cr discrepant rela¬ 
tive to the values inferred for the 2007 and 2008 lightcurves, 
while the value inferred for Rp/R^ is > 4a discrepant. The 
source of this disagreement is not clear. We also performed 
the lightcurve fitting using polynomial xy decorrelations and 


the PLD method of Deming et al. (20141, but these anal¬ 
yses gave similarly discrepant results. We therefore suspect 
that there is either an issue with the data itself or our pho¬ 
tometric reduction for this particular lightcurve. For these 
reasons, we do not consider the 2011 and 2014 3.6pm transit 
lightcurves any further in this paper. However, we consider 
our analyses for the eclipses in these lightcurves to be more 
robust, as they gave results that are consistent with those 
obtained at other epochs. 
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HD209458b transmission 






HD209458b emission 





Figure 3. Comparison of transmission Rp/Ri, (top) and emission Fp/Fi, (bottom) results obtained in the present study with those 
published in the literature. Wavelength channels are indicated in the top left corner of each axis. Independent analyses of the same 
dataset are linked by solid black lines: filled symbols show the values obtained in the current study and unfilled symbols show values 
obtained by other authors. The latter are labelled using the same abbreviations adopted in Table [T] 


For the remaining eight transit lightcurves, inferred val¬ 
ues for Rp/R^ are shown in Figure along with values 
previously published in the literature. In the 3.6^m chan- 
nel, we find Rp/R^, = 0.12077tH!5g^| for the 2007 Dec 31 
lightcurve, which is in agreement with the value of Rp/Rt, — 
0.120835 ± 0.00054 reported by Beaulieu et al. ( 20101 ). We 
obtain a somewhat higher value of Rp/R^ = 0 . 12220 ^'qqo 62 
for the 2008 Jul 19 lightcurve, compared with the value of 
Rp/Rt ~ 0.120387 ± 0.00053 obtained by Beaulieu et al. 

Our results for the 4.5 fj,m lightcurves are in good agree¬ 


ment both with each other and with the values previously 


published by Beaulieu et al. (20101 and Zellem et al. 120141. 

For the 5.8 fim 2007 Dec 31 and 2008 Jul 19 
lightcurves, we find Rp/Rt, — 0.12007lg;go265 S'lid Rp/Ri, = 
0.11880lg go 272 5 respectively. Although these values are con¬ 
sistent with each other, they are 1.7a and 1.9a lower, respec¬ 
tively, than the corresponding values of Rp/R^ = 0.1246 ± 


0.00095 and Rp/Ri, = 0.1244±0.00059 obtained by Beaulieu 
et al. (20101. In particular, our uncertainties are ~ 2.5-3 


times larger than those of Beaulieu et al. This is most likely 
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due to the flexibility of the GP models allowing broader 
ranges of function space to be marginalised over, which in 
turn maps out broader degeneracies between the transit pa¬ 
rameters and the systematics contributions. 

For the S.O/rm channel, we find good agreement between 
our inferred parameters for the 2007 Dec 24 and 2008 Jul 
22 lightcurves, with i?p/7?* = 0.12007to'ooii4 Rp/R* — 
0.1199llQ gQ 073 , respectively. For the 2008 lightcurve, our 
value is 4.7(7 lower than the value of Rp/Rt, = 0.1240 zh 


0.00046 published by Beaulieu et al. (20101. 


7.2 Emission {Fp/F^,) 

As with the transmission measurements, we find consistent 
results across epochs and wavelength channels for the in¬ 
ferred emission values (Figure [^. Our results are mostly 
in agreement with values previously published in the lit¬ 
erature, but with typically larger uncertainties. There are 
two notable exceptions. Firstly, for the 4.5/im 2005 Nov 28 
lightcurve, we obtain Fp/Fi, = 0.1451o'qi 7%> which is 2.9 (j 
lower than the value of Ep/E* = 0.213 ± 0.015% published 
by [Knutson et ah (20081, but in agreement with the value of 
Ep/E* = 0.134 ± 0.035% published by Diamond-Lowe et al. 


(20141 for the same lightcurve. Our revision brings the value 
into good agreement with those obtained in the same wave¬ 
length channel at different epochs, both in the current study 


and by Zellem et al. (20141. Secondly, for the 2005 Nov 28 
5.8/im lightcurve, we obtain Ep/E* = 0.142lQ'ggg%, which 
is 2.3(7 lower than the value of Ep/E* = 0.310±0.043% pub¬ 


lished by Knutson et al. (20081, but in agreement with the 


value of Fp/Fi, = 0.134 dz 0.035% published by 

Diamond- 

Lowe et al. 

(20141. Finally, we note that unlike 

Diamond- 

Lowe et al. 

20141, who obtained Ep/E* values for the 2005 


Nov 28 and 2007 Dec 25 8.0pm lightcurves that conflicted 
at the 3.6(7 level, our values are in good agreement with 
each other due to the relatively large uncertainty we obtain 
for the 2005 Nov 28 eclipse depth. We therefore favour the 
more conservative uncertainty estimate provided by our GP 
analysis in this case. 


7.3 Orbital parameters (a/E*, 6) 

The normalised semimajor axis a/E* and orbital inclination 
i values recovered from the primary lightcurve analyses are 
plotted in Figure]^ Note the clear correlation between both 
parameters, i.e. higher values of o/E* are associated with 
higher values of i, and vice versa. This reflects the fact that 
these parameters exert opposing influences on the transit 
duration. 

We expect a/E* and i remain constant in time across 
the different wavelength channels. Gomputing the weighted 
arithmetic mean across epochs, we obtain a/E* = 8.87 ± 
0.05, h = 0.499 ± 0.008 and i = 86.78 ± 0.07 deg. We caution 
that the quoted uncertainties for the weighted means are 
likely to be underestimated, as they are calculated by com¬ 
bining multiple measurements under the assumption that 
the errorbars are normally distributed, which is not neces¬ 
sarily true. Nonetheless, these values are consistent at the 
~ 1(7 level with the values of a/E* = 8.77 ± 0.07 and 
i = 86.76 ± O.lOdeg obtained by [Beaulieu et ’sE. (20101, 
and a/E* = 8.810lQ'Qg9 and i = 86.69/lg 'j’g deg (jbtained 
by Zellem et al. ( 2014[ ). 



Figure 4. Normalised semimajor axis a/Ri, and orbital inclina¬ 
tion values obtained from the MCMC analyses. Lightcurve epochs 
are labelled along the horizontal axis, with different marker sym¬ 
bols for each wavelength channel: purple squares for 3.6pm, green 
circles for 4.5pm, blue triangles for 5.8pm and red diamonds for 
8.0pm. Weighted mean values are labelled on each axis, and shown 
as horizontal black lines with grey shaded regions indicating the 
corresponding la uncertainties. 


Indeed, the values obtained for a/E* and b = a/R^, cost 
in the present study are overall more consistent across wave¬ 
lengths and epochs than previously published values. For 


instance, the a/E* and i values reported by Beaulieu et al. 


(20101 for the 8.0 pm channel differ from most of the values 
they obtain in the other channels by 3-5(7. The relative con¬ 
sistency of the orbit parameters derived in the current study 
therefore offers further evidence that the GP modelling ap¬ 
proach is doing a good job of accounting for the lightcurve 
systematics and providing realistic parameter uncertainties. 


7.4 Ephemeris 

Using the transit and eclipse mid-times listed in Table [^ we 
computed the ephemeris assuming a constant orbital period 
P, with the linear relation To = Tmid{n) + nP, where n is 
the number of orbital periods since the reference epoch Tq. 
Eclipse mid-times were treated as occurring precisely 0.5E 
after the immediately-preceding transit. Before performing 
the fit, we converted the BJD timestamps from Goordinated 
Universal Time (UTG) to the Barycentric Dynamical Time 


(TDB) standard, as recommended by Eastman et al. (20101. 
To do this, we added the appropriate number of leap sec¬ 
onds to the UTC timestamps, namely: 64.184 sec for the 
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Table 5. Results of the ephemeris fit to transit and eclipse mid-times assuming a constant orbital period and zero eccentricity. As 
discussed in Section |7.4| fits were performed both with and without the 2010 Jan 19 transit mid-time. As discussed in the text, we favour 
the ephemeris obtained with the 2010 Jan 19 transit excluded from the fit. 





O — C residual (min) 

Date 

Signal 

Channel (/im) 

Including 2010 transit 

Excluding 2010 transit 

2005 Nov 28 

Eclipse 

3.6 

-6.6 ±6.7 

-1.8 ±6.7 



4.5 

-3.7 ± 2.6 

1.1 ± 2.7 



5.8 

-2.4 ±8.1 

2.4 ±8.1 



8.0 

3.8 ±4.8 

8.6 ±4.8 

2007 Dec 24 

Eclipse 

8.0 

-0.5 ±0.2 

-0.1 ± 0.2 

2007 Dec 25 

Transit 

8.0 

2.4 ± 1.5 

2.8 ± 1.5 

2007 Dec 31 

Transit 

3.6 

-0.1 ± 0.4 

0.3 ±0.4 



5.8 

-0.7 ± 0.5 

-0.4 ± 0.5 

2008 Jul 19 

Transit 

3.6 

0.5 ±0.3 

-0.3 ±0.3 



5.8 

0.6 ±0.5 

-0.2 ±0.5 

2008 Jul 22 

Transit 

4.5 

1.6 ±0.4 

0.8 ±0.4 



8.0 

0.9 ±0.3 

0.1 ± 0.3 

2010 Jan 18 

Eclipse 

4.5 

5.1 ± 1.2 

1.1 ± 1.3 

2010 Jan 19 

Transit 

4.5 

-0.2 ±0.2 

-4.7 ± 0.5 

2010 Jan 21 

Eclipse 

4.5 

1.5 ±2.1 

-2.5 ± 2.1 

2011 Jan 12 

Eclipse 

3.6 

4.6 ±2.6 

-1.5 ± 2.7 

2011 Jan 16 

Eclipse 

3.6 

8.4 ± 2.9 

2.3 ±3.0 

2014 Feb 13 

Eclipse 

3.6 

16.2 ±5.8 

3.7 ± 5.9 

2014 Feb 17 

Eclipse 

3.6 

12.6 ±2.0 

0.1 ± 2.3 



Period (day) 

3.5247361 ± 6 x lO"’’ 

3.524750 ± 1 x 10"® 



To (BJDtdb) 

2454560.80567 ± 8 x 10"® 

2454560.80588 ± 8 x 10“® 


2005 lightcurves; 65.184 sec for the 2007-2008 lightcurves; 
66.184sec for the 2010-2011 lightcurves; and 67.184sec for 
the 2014 lightcurves. We report the results in Table[^ 

When all of the mid-times are included in the fit, we 
obtain To = 2454560.80567 ± 8 x 10"® BJDtdb and P = 
3.5247361 ±6 x 10“^ day, with a reduced yf = 7.0 indicating 
a poor overall fit to the data. However, we found that when 
we exclude the 2010 Jan 19 transit mid-time, we obtain To = 
2454560.80588±8 X 10"® BJDtdb and P = 3.524750± 1 x 
10“® day, with a much improved reduced yf = 1.0. It is 
not clear why our 2010 Jan 19 transit fails to fit a linear 
ephemeris, although we note that our measured mid-time 
for this lightcurve is in excellent agreement at the O.Olcr 
level with the value obtained by Zellem et al. (20141. A full 
investigation into the cause of this discrepancy is beyond the 
scope of the current paper. However, given that the other 
eighteen mid-times are well fit by a linear ephemeris, we 
conclude that the data is consistent with a constant orbital 
period. 


8 DISCUSSION 

The results outlined in Section [3 demonstrate the effective¬ 
ness of the GP modelling approach for handling systematics 
in IRAC lightcurves. Compared with those previously pub¬ 
lished in the literature, the planet properties inferred from 


the GP analyses are overall more consistent across differ¬ 
ent epochs and, for the wavelength-independent properties, 
across the different wavelength channels. In many cases, this 
is due to the GP analyses giving uncertainty estimates that 
are up to ~ 4 times larger than those reported by other 
authors. For the reasons given in Section we argue that 
the GP uncertainties provide a more realistic reflection of 
our ignorance. This is primarily because the GP models of¬ 
fer greater flexibility for handling systematics that do not 
have a well-understood functional form, compared with the 
simple parametric approximations used widely in the liter¬ 
ature. Marginalisation of the GP model posterior distribu¬ 
tions therefore allows us to more exhaustively explore possi¬ 
ble degeneracies between the planet signal and systematics, 
and incorporate these into the uncertainties associated with 
the inferred planet properties. However, recall from Section 
[^that we fixed the covariance parameters to their MLE val¬ 
ues before marginalising over the planet parameters with 
MGMC, which means the uncertainties quoted in Table 
should in fact be regarded as lower limits to the true uncer¬ 
tainties. 


Further verification of the reliability of GP models for 
inferring planet parameters from IRAG lightcurves could be 
obtained by systematically applying the method to synthetic 
datasets, similar to what was done by Gibson (20141. This 
would allow us to directly compare inferred values for the 
planet parameters with their known values. The challenge 
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with such an approach, of course, is that it requires a realis¬ 
tic simulation of the IRAC systematics, for which we do not 
possess a functional form. With this caveat in mind, how¬ 
ever, such an investigation would be useful in the future. 

8.1 Atmosphere implications 

The IRAC data analysed in the current study address 
two fundamental hypotheses concerning the nature of 
HD 209458b’s atmosphere. The first of these is the claimed 
detection of water absorption in transmission made by 
Beaulieu et al. (2010|, based on the larger values for Rp/Ri, 


that those authors measured for the 5.8/rm and S.O^m chan¬ 
nels relative to the 3.6^m and 4.5/rm channels. The second 
is the inference of a thermal inversion in the atmosphere, 
based on the measurement of deeper eclipses in the 4.5/im 
and 5.8/im channels relative to the 3.6/im channel by|Knut^ 


son et al. (20081. 

Figure shows the transmission and emission measure¬ 
ments made in the current study with overplotted model 
spectra computed using the ID radiative transfer code ATMO 
( [Amundsen et aL]|2014[ [Tremblin et al.|2015P , assuming so¬ 
lar abundances and radiative-convective equilibrium. Mod¬ 
els were generated under two different scenarios for the day- 
to-night heat redistribution efficiency: namely, uniform re¬ 
distribution and zero redistribution (respectively, cu = 0.5 
and a; = 1 in Figure |^. Number fractions as a function 
of atmospheric pressure are given for the major molecules 
in the lower left panel of Figure and the corresponding 
pressure-temperature (PT) profiles are shown in the lower 
right panel of Figure]^ 

Our transmission results are in poor agreement with 
the clear-atmosphere model predictions. In particular, we 
do not see any enhancement in the opacity at 4.5/im rela¬ 
tive to 3.6/im due to water and carbon monoxide absorp¬ 
tion. As was mentioned in Section [7.1| our transmission re¬ 
sults are also in conflict with those originally presented by 


Beaulieu et al. (20101. Specifically, Beaulieu et al. measured 


larger effective radii for the planet in the 5.8 /tm and 8.0 fim 
channels, which coincide with a water absorption band, rela¬ 
tive to the 3.6 /im and 4.5 /tm channels. We, however, obtain 
significantly larger uncertainties in the 5.8/tm channel and 
find that the effective radius is constant, or possibly decreas¬ 
ing modestly, across the 3.6-8.0/tm IRAC wavelength range. 
The difference between the results obtained in the current 
study and those of Beaulieu et al. are likely due to the two 
approaches used for treating the ramp systematics. For the 
5.8/rm channel, Beaulieu et al. truncated the first section of 
the lightcurve and fit a linear trend in time to the remainder, 
and for the 8.0 /im channel they decorrelated the ramp using 
a quadratic polynomial in logarithmic time. The GP model 
adopted in the current study should be capable of replicat¬ 
ing both these explicit functional forms, and indeed, allows 
marginalisation over an even broader function space (Section 
[^. We also confirmed that consistent results were obtained 
with the GP model when sections of varying duration were 
truncated from the start of the lightcurve. Based on our re¬ 
vised estimates of the planet’s effective radii, we conclude 
that there is no evidence for water absorption in transmis- 


of HD 209458b. Our results do not contradict this picture, 
and could suggest that the effect of the opacity inferred by 
Deming et al. at 1.4/im remains significant out to the IRAC 
wavelengths. 

We also modified the clear-atmosphere transmission 
model by fitting for an opaque cloud deck to simulate a grey 
opacity source across the IRAC wavelength range. This was 
done by allowing the clear atmosphere model to shift verti¬ 
cally while simultaneously setting the absorption to be con¬ 
stant at the cloud deck altitude. Thus our simple model had 
two tunable parameters: the overall vertical shift of the clear- 
atmosphere model and the cloud deck altitude. Our best-fit 
to the data gives a reduced of 1-4, with an opaque cloud 
deck at Rp/Ri, — 0.1210 ± 0.003. As can be seen in the top 
left panel of Figure the resulting “cloudy” model is sim¬ 
ply a horiztonal line, implying that the IRAC data show no 
evidence for even reduced-amplitude absorption features ex¬ 
tending above a cloud deck. If we instead fit a simple flat-line 
model to the data, where the only free parameter is the ver¬ 
tical level of the opaque cloud deck, the number of degrees of 
freedom increases from six to seven, and the reduced im¬ 
proves to 1.2. Due to the size of our uncertainties, however, 
we cannot rule out muted absorption features of similar am¬ 


plitude to the H 2 O feature measured by Deming et al. (2013 \ 


at 1.4/tm with WFC3. Furthermore, we note that our result 
is consistent at the ~ la level with the cloud deck altitude 
implied by the WFC3 transmission spectrum. Such a com¬ 
parison should be treated with caution, however, as Deming 
et al. fixed a/R* = 8.95 and i = 86.93° in their lightcurve 
fits (taken from [Knutson et al.||2007b l, which would intro¬ 
duce an offset in the absolute level those authors derive for 
Rp/Ri, relative to our study. 

For the emission, we produced models both with and 
without inverted pressure-temperature profiles (lower right 
panel of Figure]^. To produce thermal inversions in the for¬ 
mer models, we artificially included titanium oxide (TiO) 
and vanadium oxide (VO) with a constant abundance 
throughout the atmosphere (lower le ft panel of Figure [^. 
Our emission results reinforce those of IDiamond-Lowe et al.l 


(20141, who argued that there is no evidence for a thermal 


inversion in the atmosphere of HD 209458b based on the 
revised estimates for the 4.5/im and 5.8/im eclipse depths. 
This can be seen in Figurej^ where the models with a ther¬ 
mal inversion are shown to provide a very poor match to 
the data. For the models without a thermal inversion, the 
match is better, but still poor. In particular, the latter mod¬ 
els underpredict the emission in the 4.5/im channel, where 
an absorption feature due to carbon monoxide is expected 
to block radiation emitted from the planetary atmosphere. 
The fact that we do not detect this absorption is in line 
with the recent non-detection of carbon monoxide in the 


sion over the IRAG bandpasses. Indeed, Deming et al. ( 2013[ ) 
have measured a muted water feature at 1.4/im using WFC3, 
and interpreted this as evidence for haze in the atmosphere 


dayside hemisphere made by Schwarz et al. (20151 using 
high-resolution spectroscopy. 

We also fit the emission data with a model that assumes 
the planet radiates as an isothermal blackbody. To gener¬ 
ate this model, an ATLAS spectrum ( [Kurucz I pT^ pMt 
computed specifically for HD 209458 was used for the stellar 
emissiorj^ and the radius ratio was fixed to Rp/R* = 0.121. 
Assuming a Planck spectrum for the planet, a temperature 


® http://kurucz.harvard.edu/stars/hd209458 
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Number fraction 




Temperature 


Figure 5. Top panels show the transmission Rp/Ri, (left) and emission Fp/F* (right) measurements obtained in the current study for 
HD 209458b, with colours indicating wavelength channels. IRAC bandpasses are indicated at the bottom of both axes as solid black lines. 
For measurements obtained in the same wavelength channel, small horizontal offsets have been applied for clarity. Both panels show 
ID ATMO models assuming solar abundances and chemical equilibrium. Corresponding chemical abundances are shown in the lower left 
panel and pressure-temperature profiles are shown in the lower right panel. For the abundances, solid lines correspond to uniform heat 
redistribution from the dayside to the nightside (i.e. ui = 0.5), while dashed lines correspond to zero heat redistribution (i.e. ui = 1). For 
the transmission, the light blue line shows the clear atmosphere obtained assuming uj = 0.5, while the grey line shows a simple model 
corresponding to an opaque cloud deck at i?p/i?* = 0.1210, based on fitting a simple cloud model to the data as described in Section 
|8.1| For the emission, clear atmosphere models are shown for uj = 0.5 and ui = 1, both with and without TiO/VO in the atmosphere. For 
those models including TiO/VO, the species were assumed to be mixed uniformly throughout the atmosphere, as shown by the brown 
line in the lower left panel. The inclusion of TiO/VO results in inverted temperature-pressure profiles, as shown in the lower right panel. 
Also shown in the upper right panel is the emission spectrum obtained using an ATLAS stellar model for HD 209458 and assuming the 
planet radiates as a blackbody with a temperature of 1484 K. In the top panels, square symbols give the model values integrated over 
the corresponding IRAC bandpasses. 


of 1484 ± 18K was found to giv e the best fit to the data, 

(2014 1 have 


Hansen et al. 


with a reduced x of 1-5- Indeed, 
recently put forward the case that blackbody radiation can 
explain the majority of IRAC emission data that has been 
pnblished for exoplanets to date, largely dne to underesti¬ 
mated uncertainties for the eclipse depths. The acceptable fit 
provided by the isothermal blackbody for onr data supports 
this hypothesis in the case of HD 209458b, especially consid¬ 


ering the reduced could be even lower if we marginalised 
over the covariance parameters and obtained larger nncer- 
tainties (Section]^. 

One possibility is that the emission measurements are 
probing an isothermal layer of the atmosphere above a clond 
deck that extends across the dayside hemisphere. This could 
simultaneously explain the lack of absorption features de¬ 
tected in the IRAC transmission data, and the muted H 2 O 
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feature detected in the WFC3 transmission data, discussed 
above. We stress, however, that this is a speculative scenario, 
poorly constrained by the existing observations. 


9 CONCLUSION 

We have presented an analysis of IRAC transit and eclipse 
lightcurves for HD 209458b. By binning the lightcurves in 
time, it was possible to perform the lightcurve analyses 
using GP models. GPs allow transit lightcurve models to 
be elegantly defined within a rigorous Bayesian framework. 
They provide a natural mechanism for handling poorly- 
understood systematics in the data that are unrelated to the 
astrophysical signal of interest. Uncertainty is propagated 
through all levels of the model in a clear and transparent 
manner, and Occam’s razor is automatically implemented, 
mitigating against overfitting. 

We have made a number of significant revisions to pre¬ 
viously published results, and in many cases the uncertain¬ 
ties for inferred planet parameters have been increased by 
factors of ~l-4. This can largely be attributed to the flex¬ 
ibility of the GP models, which allow complex correlations 
to be handled with a small number of free parameters. The 
latter point is important, as it means that marginalisation 
over the model parameter space remains computationally 
tractable, allowing uncertainties that realistically quantify 
the degeneracies between the planet signal and instrumen¬ 
tal systematics to be derived. 

We obtain an overall improvement in consistency for the 
normalised semimajor axes a/7?* and orbital inclinations i 
across different epochs and wavelength channels, compared 
to results that have been previously published for these 
datasets. This provides evidence that the GP models are 
effectively accounting for the systematics, and typically less 
prone to underestimating uncertainties compared with other 
lightcurve fitting approaches used in the literature. 

The revised GP analyses presented here draw into ques¬ 
tion a number of claims that have previously been made re¬ 
garding the atmosphere of HD 209458b, including the detec¬ 
tion of water absorption in transmission and the inference 
of an inverted PT profile for the dayside hemisphere. In¬ 
stead, our transmission measurements are consistent with a 
featureless spectrum, and our emission measurements are fit 
reasonably well assuming the planet radiates as an isother¬ 
mal blackbody with a temperature of 1484 ± 18 K. 

Taken together, our results illustrate how sensitive 
IRAC analyses are to the systematics treatment. GP anal¬ 
yses have been shown to produce results that are gener¬ 
ally more stable, and with uncertainties that are relatively 
conservative, compared to those obtained using other com¬ 
mon approaches. However, we do tend to find good agree¬ 
ment with results obtained using pixel-mapping techniques 
in the 3.6/im and 4.5/rm channels, which is unsurprising 
given that GPs are in essence quite similar to pixel-mapping. 
Nonetheless, many results have been published using simpler 
xy polynomial decorrelations and parametric ramp models, 
which may not always be adequate. Combined with the lack 
of spectral resolution afforded by the broad bandpasses, our 
results suggest that statements made previously in the lit¬ 
erature about exoplanet atmospheres relying heavily on the 


interpretation of IRAC data should be regarded with cau¬ 
tion. 
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Table Al. Maximum likelihood estimates for the GP covariance 
parameters. These are the values that were fixed for the MCMC 
analyses described in Section]^ 


Channel 

Date 

Signal 

■^xy 

Lx 

Ly 

At 

Lt 


(pm) 



(%) 

(pix) 

(pix) 

(%) 

(min) 

(ppm) 

3.6 

2005 Nov 28 

Eclipse 

2.0552 

0.642 

0.528 

0.0444 

37.117 

1541 


2007 Dec 31 

Transit 

0.4874 

0.169 

0.187 

0.0216 

20.708 

2185 


2008 Jul 19 

Transit 

0.3089 

0.138 

0.136 

0.0677 

0.319 

2066 


2011 Jan 12 

Eclipse (1st) 

1.0988 

0.176 

0.091 

0.0261 

5.471 

499 


2011 Jan 14 

Transit 

0.6342 

0.127 

0.145 

0.1754 

0.322 

344 


2011 Jan 16 

Eclipse (2nd) 

1.8569 

0.235 

0.233 

0.0331 

6.710 

541 


2014 Feb 13 

Eclipse (1st) 

1.5366 

0.167 

0.135 

0.0398 

5.574 

553 


2014 Feb 15 

Transit 

1.6826 

0.195 

0.196 

0.0488 

0.553 

660 


2014 Feb 17 

Eclipse (2nd) 

1.2590 

0.151 

0.118 

0.0234 

3.089 

453 

4.5 

2005 Nov 28 

Eclipse 

0.8656 

0.238 

0.254 

0.8796 

2386.526 

1391 


2008 Jul 22 

Transit 

0.0001 

143.003 

109.329 

0.0302 

5.887 

3148 


2010 Jan 18 

Eclipse (1st) 

0.2467 

0.067 

0.090 

0.0211 

25.212 

578 


2010 Jan 19 

Transit 

0.8769 

0.255 

0.231 

0.0086 

24.042 

569 


2010 Jan 21 

Eclipse (2nd) 

0.8468 

0.235 

0.270 

0.0182 

15.931 

566 

Channel 

Date 

Signal 

A, 

Lr 

h 

At 

Lt 


(pm) 



(%) 

(log [min]) 

(min) 

(%) 

(min) 

(ppm) 

5.8 

2005 Nov 28 

Eclipse 

0.0732 

6.187 

2796.4 

0.0575 

195.8 

3956 


2007 Dec 31 

Transit 

0.3080 

6.127 

489.8 

0.0000 

67.9 

3037 


2008 Jul 19 

Transit 

0.7993 

12.986 

0.1 

0.0000 

> 10^ 

2986 

8.0 

2005 Nov 28 

Eclipse 

0.2781 

6.277 

1142.0 

0.0000 

> lO"' 

3749 


2007 Dec 24 

Transit 

0.0000 

357.352 

291853.4 

0.2427 

1221.2 

1163 


2007 Dec 25 

Eclipse 

0.0163 

5.887 

12729.4 

0.0000 

> 10^ 

1215 


2008 Jul 22 

Transit 

1.8164 

9.791 

1.0 

0.0000 

> 10^ 

2213 


APPENDIX A: COVARIANCE PARAMETER 
MAXIMUM LIKELIHOOD ESTIMATES 

Table [XT] reports the maximum likelihood estimates for the 
GP covariance parameters that were fixed for the MCMC 
analyses, as described in Section 

This paper has been typeset from a fOiiX/ file prepared 

by the author. 


© 2015 RAS, MNRAS 000,[iHT8l 







